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We study the spectral width as a function of the external momentum for the dynamical struc- 
ture factor of a disordered harmonic solid, considered as a toy model for supercooled liquids and 
glasses. Both in the context of single-link coherent potential approximation and of a single-defect 
approximation, two different regimes are clearly identified: if the density of states at zero energy is 
zero, the Rayleigh law is recovered for small momentum. On the contrary, if the disorder induces 
a non vanishing density of states at zero energy, a linear behaviour is obtained. The dynamical 



structure factor is numerically calculated in lattices as large as 96"^, and satisfactorily agrees with 



a^ : 

ON . 

^ the analytical computations. 

^ . PACS numbers: 61.43.Fs,63.50.+x 

> ; 

, O I The spectrum of vibrational excitations of supercooled liquids and glasses is attracting a great deal of attention from 
' the experimental side |Q] , the numerical simulations |^ and also from the analytical point of view . Vibrational 
Q.^ I excitations on the GHz range develop in the supercooled liquid, by the same temperatures where anomalous behaviour 
. of the specific-heat is found. Some rather substance-independent features in this vibrational spectrum are found. For 
example, the vibrational density of states presents an excess respect to the usual Debye behaviour {g{uj) oc w^) known 
^-H as the Boson Peak. The dynamical structure-factor, S'(p, w), reveals well-defined sound-like peaks for wave-lengths 
^ not much greater than the inter-particle distance. Moreover, the speed of this high-frequency sound is close to the 
one of the low- frequency typical one. For fixed external momentum, p, the width of the spectralpeak grows as p , 
which has recently been recovered in the mode coupling approximation (see [Q). However, in Ref. g, a simpler model 
of a disordered three-dimensional harmonic solid was studied. There it was claimed that if some spring-constants are 
. allowed to have a small negative value, (but not so much that negative-energy modes appear) a boson peak develops. 
^\ ' In the context of the CPA approximation ||], the usual width of S{p, w) arisen from Rayleigh-scattering is reported. 
On ' It was, however, noticed that a rather different scaling appeared at the characteristic frequencies of the Boson-Peak. 

Furthermore, in Ref. |Q, it was shown a broadening for a one-dimensional disordered harmonic solid model. In 
^ , this paper, we want to investigate the problem, both by numerical and analytical means. Our main finding will be 
^ ■ that the Rayleigh-scattering p'^ broadening holds, unless the solid becomes unstable. That is, if the dynamical-matrix, 
to be defined later, has an extensive number of negative eigenvalues, the system do present sound-like peaks in its 
S{p, oj), but with a width proportional to p, at least for very small values of p. Otherwise, the standard p"^ behaviour 
Q is to be expected. 
^ . To be more specific, our dynamical matrix in D dimensions is given by 



H =^^(t)x'Hxy(f>y , (1) 



xy 

D 



+ ^^^f^{Sxy-Sx,y-,), (2) 

where fl is the lattice unit-vector in the /i direction. Notice that we cannot separate the longitudinal and transverse 
modes since, for the sake of simplicity, our matrix has no internal indices as in Ref. The a^.^ is the random part 
of the spring constant that joins the sites y and y + p.. In a finite lattice, periodic boundary conditions are applied. 
The dynamical matrix, verifies the constraint related with traslational symmetry: the vector of all equal components 
is eigenvector with zero eigenvalue. For a harmonic solid, the object that we will study is the dynamical structure 
function, S{p, E), in the energy rather than in the frequency-domain: 



S{p,E)^~-Im\im{p\ \p). (3) 

In the above expression, the bra-ket notation has been used and \p) stands for a normalized plane-wave of momentum 
p. As usual, the overline represents the average over the random variables a. We choose to work in the energy, rather 
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than in the frequency domain, since we mainly want to consider the case where a significant fraction of the spectrum 
is negative, and so the transformation E = uj'^ is no longer well defined. 

We study the case where the aj,^^ are uncorrelated, random variables whose probability distribution is 

pia)^{l-p)dia)+pfia),0<p<l. (4) 

In the above equation, p is the probability of finding one defect, while / is a continuous probability function, that 
we take flat between A and (A < 0). Notice that the disorder of Ref. ||] is recovered, by taking p = 1 and / 
Gaussian. The rationale for choosing this kind of disorder is that a spring of an unusually large negative spring- 
constant induces a negative energy eigenvalue of the dynamical matrix Ti,. In fact, if one keeps in Eq.(^) only the 
spring connecting sites 1 and 2, corresponding to the large negative defect, the eigenvectors are easily shown to be 
(f'x = {Sx,i + Sx.2) I \/2 and (8x.\ — ^x.i) I V2. However, assuming that the surrounding springs have their ordered value, 
one finds that the Hamiltonian (|^) is negative for the displacement configuration (J^, ! — 5x^2)1 if the spring-constant 
verifies a < — (£> + 5). Therefore, one has a contribution of order p to the density of states over the negative spectrum. 
Moreover, if we assume that some of the surrounding spring have very small positive values, a contribution of higher- 
order in p to the density of states is generated, no matter how small is the negative value of the spring constant, 
a + 1. Therefore, if the probability of a < —1.0 is non- vanishing, the hybridization of the localized field-configuration 
described above, with the plane-waves eigenstates will non-trivially modify the eigenvectors of Ti. To study this, we 
shall perform a single-defect calculation of the resolvent. 

In this way, we learn that the order p threshold for the presence of negative eigenvalues is not a = —D — 1/2, as 
roughly shown in the introduction, but —D. This threshold separate two well defined scaling limits for the width of 
the S{q,E) at small q: 

• If the density of states is null at zero energy the imaginary part of the self-energy E is proportional to 
(Rayleigh scattering, that yields a spectral width cx p^'^^ in the frequency domain). 

• When an extensive number of negative eigenvalues is present, S cx (or p for the width in the frequency domain). 
Of course, the p^~^^ contribution will be still present but sub-dominant at low momentum. A crossover might 
be visible, depending on the strength of the disorder. 

The dynamical matrix in the presence of an unique defect of amplitude a can be split in two terms: 

Ti-xy — Tixy + TZxy (5) 

where the pure crystal matrix TiP has plane waves eigenvectors with eigenvalues 

Eo{p) = Y.^l-cosp,). (6) 

1/ 

The perturbation, hence, connects the two sites and + u: 

TZxy^ay"-'-'\y";ti){y";p\, (7) 

0_ \ ^x,y^ ^x,y^-\-i^ 



The propagator can be written as 



\v»U')s (8) 



z-H z-W z-W z-W 
where the resummation of the harmonic series gives 

T^|y°;^)-^-^(y°;H, (10) 
1 — aa(z) 

aiz) — — / - — -— . (11 

Notice that the correction term in Eq. (p^) , has an isolated singularity for the value of z satisfying 1 — aa{z) = 0. 
This value decrease monotonically from z = 0, for a = —D, and will be called Za- This isolated singularity, correspond 



with 
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to the ground state of the dynamical matrix, the residue being ^p{x)'^*p{y). Therefore, one obtain the wave- function 
for this eigenvalue 

(X / — ^ -— e^-^- , (12) 

J {2tt)^ Za - Eo[q) 

that has a localization length of order \za\~^^^. We thus see, that unless a was exceedingly close to the critical 
value —D, the eigenvector is strongly localized around the defect. The single-defect approximation to the self-energy, 
Eg. (p^ , amounts to consider that each defect only contributes to its own localized eigenvector, and that no other 
defect is within its localization length. Let us turn back to the original problem, with an extensive number of defects, 
and repeat the above calculation, neglecting all terms which contain two different defects. The matrix T is now a 
sum of terms like the one in Eq.(|o|). If we now perform the average over the a's and apply Dyson resummation, the 
self-energy part of the propagator turns out to be 



E{z,p) = pEoip) / dafia)- — , (13) 

J 1 - aa{z) 

at first order in p (the defect interactions will generate the order and higher order corrections to the single defect 
result). 

The width of the S{q, E) is simply given by the value of the imaginary part of the self-energy at the peak, whose 
position can be obtained from the real part of the self-energy, i?™'*'^(p) k, Eo{p) + ReI](Eo(p), p). Notice that our self 
energy is proportional to Eq{p) so that we are basically getting a finite renormalization of the speed of sound. In order 
to estimate the imaginary part of the self energy, it is useful to realize that for small positive values of the energy, one 
has 

'^(^ + ^~15- ' 20/2j/2-iriD/2) ■ ('^^ 

Therefore, if the probability density /(a) do not allow a to be smaller than —D, the only imaginary term in Eq.(p^) 
comes from a{E + ie), and it is of order p^E^/^ , yielding a value of order p^ at the peak. On the other hand, if a can 
be smaller than —D, an imaginary part of order p^f{—D) arises from the pole. 

Led by the functional form of the self-energy in Eq.(p^, one can consider self-energies of the form /(i?)i?o(p), also 
when p is not small, and the single defect approximation no longer holds. This is the idea lying under the well-known 
CPA approximation M, where one sets Ti{z,p) — (T{z) — l)Eo{p). A self consistency equation can be readily written: 



T{z)-a{z/T{z)){l + a^T{z)) 



-0. (15) 



It is clear that the width of the S{q,E) critically depends on the value of r(0 + ie). For the flat distribution of a 
introduced in Eq.(0), one can solve Eq.(|l5|) in the hmit of small p, as r(0) ~ 1 + bp + 0{p^) obtaining 



h = D[l-D\og + tT:^6{D + A)^ , (16) 

where 9{x) = 1 for a; < and zero otherwise. It is easy to check that the single-defect result is exactly the same as 
in Eq.(|l|). 

Fixing from now on Z? = 3, Eq.(^5|) in the general case, can be numerically solved. We choose to write it as a 
fixed-point equation, and solve it recursively. The only tricky part, is the evaluation of a{z) defined in Eq.(|ll]), that 
we make by calculating the unperturbed density of states by a Monte Carlo simulation. In this way, we are able to 
obtain estimates of T{z) with an accuracy of 10^'*. We find that in the CPA approximation, there are the same two 
regimes as in the single-defect computation, separated by a critical line, whose coordinates are shown in table |. The 
appearance of the p^ regime coincides, and it is due to, the vanishing of the negative-energy spectrum. However, as 
explained in the introduction, one expects rather that what vanish is the order p contribution to the density of states 
in this region. In fact, we expect a non- vanishing density of states all the way down to A — —1, where no negative 
spring exists. It is remarkable that when p = 1, the CPA approximation keeps a non-vanishing fraction of negative 
eigenvalues up to A = —1.2, quite close to the correct value. 

The agreement between the CPA approximation and the results from numerical simulations for the S{E,p) turns 
out to be better than 5% in the two extreme cases (see figs. |[^). That is the one where there is no negative springs 
(A = —1) and the one where the spring constant can be very negative (A = —10). Close to the CPA critical line (in 
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the unstable crystal side) the agreement is still quite good. On the (CPA) stable border side, a scaling between 
and is found, depending on the density of defects, p. For example, with the value of the density of defects p = 0.1, 
in fig. 1^) the transition between the two regimes found by the CPA approximation at A = —2.15 is very evident. 
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FIG. 1. The spectral width E (see Eq.(|l7|)) as a function of the external momentum in the CPA approximation, for a density 
of defects p = 0.1 

We numerically computed the structure factor S{E,q) for our model utilizing the method of moments jsj, which 
allows to study the statistical properties of the eigenvalues and the eigenvectors of large dynamical matrices, avoiding 
their diagonalization. This method is a clever modification of the Lanczos method and it shares the same weakness, 
namely the lack of orthogonality when a too large number of moments is computed. Another limitation is the necessity 
of setting a finite value of e in Eq.(||). A reasonable value for e is the mean distance, between eigenstates, which is 
roughly given by {2D — X)/V, where V is the lattice volume (that is 10~^ in our 96^ lattice). If the width of the 
peak is comparable with e, the results will be definitely affected by finite-size effects. The limitation related with the 
number of moments, is not serious for the central part of the S{q,E) curve, but can be rather strong if one wants 
to calculate the tails of the distribution. In practice, we have used 30 moments for 10 different disorder realizations, 
finding very satisfying results, unless the peak height approaches values of order 10^, when finite size effects turns out 
to be important. We fit our results to the Breit-Wigner form: 



SiE,q) 



(17) 



which satisfactorily describes the peak in all cases, although usually overestimates the tails of the distribution. The 
position of the peak Eg is linear in p'^ as expected. More importantly, the results from numerical simulations (fig. 
^ show very clearly that the critical line previously discussed is an artifact of the CPA approximation and a real 
system actually becomes unstable when there is an extensive number of negative springs, no matter how small is the 
fraction of them, as expected from the simple analytic considerations sketched previously. For larger values of p the 
p^ contribution is no longer negligible and a complicated intermediate scaling appears. 



4 




5 



1000 




FIG. 3. The S{q, E) function, for p = 2%/37rn/96 for n = 1, 2, 6, 10, 14, from left to right, both in the CPA approximation 
and in the numerical calculation, for p = 0.1 and A = —1.0. 
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FIG. 4. E as a function of the external momentum in the numerical simulations, for a density of defects p — 0.1. The errors 
are smaller than the data points. 

In this work, we used the fact that both the single-defect approximation and the CPA approximation |^ allow to 
write the self-energy of a disordered harmonic solid, for small E and p, as f{E)Ea{p). An unavoidable consequence 
of this functional form is that violations of the Rayleigh scaling of the width of the peaks of the S{p, E) are 
intrinsically tied with the appearance of a negative energy spectrum (z. e. solid instability). We have argued that an 
extensive number of negative eigenvalues should appear, as soon as the spring-constants are allowed to be negative. 
Our numerical calculations on 96^ disordered lattices confirm the above picture. For wave-lengths in the range [a, 10a] 
(a being the lattice constant), a crossover regime may come out, depending on the strength of the disorder, from the 
competition between the term and p^ term in the self-energy. Therefore, the vibrational excitations of the disordered 
solid can look qualitatively similar to the finite temperature Instantaneous Normal Modes of supercooled liquids and 
glasses, where negative eigenvalues are always present. 

We gratefully acknowledge interesting discussions with A. Gonzalez, M. Mezard, G. Ruocco and G. Viliani. V.M.M. 
is a M.E.C. fellow and has been partially suported by CYCyT(AEN97-1708 and AEN99-1693). Our numerical 
computations have been carried out on the Kalix2 pentium cluster of the University of Cagliari. 



TABLE I. The critical line \c{p) in the CPA approximation. 



p 


0.0001 


0.1 


0.2 


0.4 


0.6 


0.8 


1.0 


-A 


2.995(5) 


2.15(5) 


1.84(1) 


1.55(2) 


1.38(1) 


1.28(1) 


1.22(1) 
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